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Abstract 

This paper concerns the harmonic shift-invert residual Arnoldi (HSIRA) and Jacobi- 
Davidson (HJD) methods as well as their refined variants RHSIRA and RHJD for the 
interior eigenvalue problem. Each method needs to solve an inner linear system to expand 
the subspace successively. When the linear systems are solved only approximately, we are 
led to the inexact methods. We prove that the inexact HSIRA, RHSIRA, HJD and RHJD 
methods mimic their exact counterparts well when the inner linear systems are solved 
with only low or modest accuracy. We show that (i) the exact HSIRA and HJD expand 
subspaces better than the exact SIRA and JD and (ii) the exact RHSIRA and RHJD 
expand subspaces better than the exact HSIRA and HJD. Based on the theory, we design 
stopping criteria for inner solves. To be practical, we present restarted HSIRA, HJD, 
RHSIRA and RHJD algorithms. Numerical results demonstrate that these algorithms are 
much more efficient than the restarted standard SIRA and JD algorithms and furthermore 
the refined harmonic algorithms outperform the harmonic ones very substantially. 

Keywords. Subspace expansion, expansion vector, inexact, low or modest accuracy, 
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1 Introduction 

Consider the linear eigenproblem 

Ax = Xx, x H x = 1, (1) 

where A G c nxn is large and possibly sparse with the eigenvalues labeled as 

< |Ai — a\ < | A2 — er| < • • • < |A n — a\, 

where a G C is a given target inside the spectrum of A. We are interested in the eigenvalue 
Ai closest to the target a and/or the associated eigenvector x\. This is called the interior 
eigenvalue problem. We denote (Ai,xi) by (A, x) for simplicity. 

The interior eigenvalue problem arises from many applications pQ. Projection methods 
have been widely used for it [2, 19,20,25,26). For a given subspace V, the standard projection 
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method, i.e., the Rayleigh-Ritz method, seeks the approximate eigenpairs, i.e., the Ritz pairs, 
(u, y) satisfying 

(A - ul)y IV, y € V. (2) 

As is well known, the standard projection method is quite effective for computing exterior 
eigenpairs, but is not for computing interior ones. For interior eigenpairs, one usually faces 
two difficulties when using the standard projection method: (i) the Ritz vectors might be 
very poor [HE] even if the subspace that contains good enough information on the desired 
eigenvectors; (ii) even though there are Ritz values which are good approximations to the 
desired eigenvalues, it is usually hard to recognize them and to select the correct Ritz pairs 
to approximate the desired interior eigenpairs [171118] . Because of possible mismatches, the 
Ritz pairs may misconverge or converge irregularly as subspace V is expanded or restarted. 

To overcome the above second difficulty, one popular approach is to apply the standard 
projection method to the shift-invert matrix (A — a/) -1 , which transform the desired eigen- 
values of A into the dominant (exterior) ones. However, if it is very or too costly to factor 
A — al, as is often the case in practice, one must resort iterative solvers to approximately 
solve the linear systems involving A — al. Over years, one has focused the inexact shift- 
invert Arnold (SIA) type methods where inner linear systems are solved approximately. It 
has turned out that inner linear systems must be solved with very high accuracy when ap- 
proximate eigenpairs are of poor accuracy [5l[2Tl[22]. As a result, it may be very expensive 
and even impractical to implement SIA type methods. 

As an alternative, one uses the harmonic projection method, i.e., the harmonic Rayleigh- 
Ritz method, to seek the approximate eigenpairs, the harmonic Ritz pairs, (z/, y) satisfying 

(A - (ul)y ±(A- aI)V, y E V. (3) 

The harmonic projection method has been used widely for the interior eigenvalue problem 
[2j[20l[25l[26]. The method can be derived by using the standard Rayleigh-Ritz method on 
(^4 — al)" 1 with respect to the subspace (^4 — aI)V without factoring A — al. The harmonic 
projection method has the advantage that is the Ritz value of (^4 — al)' 1 . This is 
helpful to select the correct approximate eigenvalues since the desired interior eigenvalues 
near a have been transformed into the exterior ones that are much easier to match. So the 
harmonic projection method overcomes the second difficulty well. 

However, the harmonic projection method inherits the first difficulty since the harmonic 
Ritz vectors may still be poor, converge irregularly and even fail to converge; see [12] and 
also [25^126] for a systematical account. To this end, the refined projection method, i.e., 
the refined Rayleigh-Ritz method, has been proposed in [3115] that cures the possible non- 
convergence of the standard and harmonic Ritz vectors. The refined method replaces the 
Ritz vectors by the refined Ritz vectors that minimize residuals formed with the available 
approximate eigenvalues over the same projection subspace. The refined projection principle 
has also been combined with the harmonic projection method. Particularly, the harmonic and 
refined harmonic Arnoldi methods were proposed in [10] that are shown to be more effective 
for computing interior eigenpairs than the standard and refined Arnoldi methods. Due to the 
optimality, refined (harmonic) Ritz vectors are better approximations than (harmonic) Ritz 
vectors. The convergence of the refined projection methods has been established in a general 
setting [ISICES]; see also [251126]. 

A necessary condition for the convergence of any kind of projection methods is that the 
subspace contains enough accurate approximations to the desired eigenvectors. Therefore, a 
main ingredient for the success of a projection method is to construct a sequence of subspaces 
that contain increasingly accurate approximations to the desired eigenvectors. To achieve this 
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goal, it is naturally hoped that each expansion vector makes contribution to the desired 
eigenvectors as much as possible. It turns out that preconditioning plays a key role in 
obtaining good expansion vectors [18l[23]. Such kind of methods includes the Davidson 
method [HPI], the Jacobi-Davidson (JD) method [23] and shift-invert residual Arnoldi (SIRA) 
method [16]. 

The Davidson method can be seen as a preconditioning Arnoldi (Lanczos) method and a 
prototype of the SIRA and JD methods, in which an inner linear system involving A — al 
needs to be solved at each outer iteration. When a factorization of A — al is impractical, only 
iterative solvers are viable for solving the inner linear systems. This leads to the inexact SIRA 
and JD. A fundamental issue on them is how accurately we should solve inner linear systems 
in order to make the inexact methods mimic the exact counterparts well, that is, the inexact 
and exact ones use comparable outer iterations to achieve the convergence. For the simplified 
or single vector JD method without subspace acceleration, it has been argued that it may not 
be necessary to solve inner linear system with high accuracy and a low or modest accuracy 
may suffice [6ll23j. This issue is more difficult to handle for the standard JD method. Let the 
relative error of approximate solution of an inner linear system be e. Lee and Stewart [16j 
have given some analysis of the inexact SIRA under a series of assumptions. However, their 
main results appear hard to justify or interpret, and furthermore no quantitative estimates 
for e have been given eventually. Using a new analysis approach, the authors |14| have 
established a rigorous general convergence theory of the inexact SIRA and JD in a unified 
way and derived quantitative and explicit estimates for e. The results prove that a low or 
moderate e, say [10 -4 , 10~ 3 ], is generally enough to make the inexact SIRA and JD behave 
very like the exact methods. As a result, the methods are expected to be more practical than 
the inexact SIA method. The authors have confirmed the theory numerically in [14]. 

Because of the merits of harmonic and refined harmonic projection methods for the interior 
eigenvalue problem, in this paper we consider the harmonic and refined harmonic variants, 
HSIRA, RHSIRA and HJD, RHJD, of the standard SIRA and JD methods. Exploiting the 
analysis approach and some results in [14], we establish the general convergence theory and 
derive the estimates for e for the four inexact methods. The results are similar to those for 
the SIRA and JD methods in [14] . They prove that, in order to make the inexact methods 
behave like their exact counterparts, it is generally enough to solve all the inner linear systems 
with only low or modest accuracy. Furthermore, we show that the exact HSIRA and HJD 
expand subspaces better than the exact SIRA and JD and, in turn, the exact RHSIRA 
and RHJD expand subspaces better than the exact HSIRA and HJD. These results and 
the merits of the harmonic and refined harmonic methods mean that the harmonic Ritz 
vectors and refined harmonic Ritz vectors are better for subspace expansions and restarting 
than the standard Ritz vectors. Based on the theory, we design stopping criteria for inner 
solves. To be practical, we present restarted HSIRA, RHSIRA, HJD and RHJD algorithms. 
We make numerical experiments to confirm our theory and demonstrate that the restarted 
HSIRA, HJD and RHSIRA and RHJD algorithms are considerably more efficient than the 
restarted standard SIRA and JD algorithms and furthermore the restarted RHSIRA and 
RHJD outperform the others very greatly for the interior eigenvalue problem. 

This paper is organized as follows. In Section [21 we describe the harmonic and refined 
harmonic SIRA and JD methods. In Section [3J we present our results on e. In Section |4] we 
describe restarted algorithms. In Section [5l we report numerical experiments to confirm our 
theory. Finally, we conclude the paper and point out future work in Section [6] 

Some notations to be used are introduced. Throughout the paper, denote by || • || the 
Euclidean norm of a vector and the spectral norm of a matrix, by / the identity matrix with 
the order clear from the context, and by the superscript H the complex conjugate transpose 
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of a vector or matrix. We measure the deviation of a nonzero vector y from a subspace V by 

\\y\\ 

where Py is the orthogonal projector onto V. 

2 The harmonic SIRA and JD methods and their refined ver- 
sions 

Let the columns of V form an orthonormal basis of a given general subspace V and define 
the matrices 

H = V H {A-aI) H V, (4) 
G = V H (A-aI) H (A-aI)V. (5) 

Then the standard projection method computes the Ritz pairs (v, y) satisfying 

H z = Liz with ||z|| = 1, 

u = fi + a, (6) 

y = Vz, 

and selects z/s closest to cr, which correspond to the smallest //'s in magnitude, and the 
associated y's to approximate the desired eigenpairs of A. Given an initial subspace V and 
the target a, we describe the standard SIRA and JD methods as Algorithm 1 for our later 
use. 

Algorithm 1 The SIRA/Jacobi-Davidson method with the target a 

Given the target a and a user-prescribed convergence tolerance tol, suppose the columns 

of V form an orthonormal basis of the initial subspace V. 

repeat 

1. Use the standard projection method to compute the Ritz pair (y,y) of A with 
respect to V, where v = A and \\y\\ = 1. 

2. Compute the residual r = Ay — vy. 

3. In SIRA, solve the linear system 

(A - al)u s = r. (7) 
In JD, solve the correction equation 

{I -yy H ){A-aI){I -yy H )uj = -r (8) 

for uj A. y. 

4. Orthonormalize the us or uj against V to get the expansion vector v. 

5. Subspace expansion: V = [ V v J . 
until llrll < tol 



The harmonic projection method seeks the harmonic Ritz pairs (fx, y) satisfying ([3]), which 
amounts to solving the generalized eigenvalue problem of the pencil (H, G): 



Hz 


= ±Gz 






V 


= li + a 


y 


= Vz, 
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and selects z/s closest to a, which corresponds to the smallest p's in magnitude, and the 
associated y's to approximate the desired eigenpairs of A. 

Alternatively, we may use the Rayleigh quotient p of A with respect to the harmonic Ritz 
vector y to approximate A as well. Recall the definition ([4]) of H and note from ([9]) that 
y = Vz and \\z\\ = 1. Then 



y Ay _ z»V H (A-aI + aI)Vz _ HrjH 

y H y 



P=-7M7r = z h v h Vz =z H H«z + a, (10) 



so we can compute p very cheaply. It is proved in [12] that once a is very close to the desired 
eigenvalue A, then the harmonic projection method may miss A and the harmonic Ritz value 
fj, may misconverge to some other eigenvalue; see \12\ Theorem 4.1]. However, whenever y 
converges to the desired eigenvector x, the Rayleigh quotient p must converge to A no matter 
how close a is to A; see [HI Theorem 4.2]. Therefore, it is better and safer to use p, rather 
than p, as an approximate eigenvalue. Another merit is that p is optimal in the sense of 

p = arg min \\(A — ul)y\\. 



In the sequel, as was done in the literature, e.g., [12|[T8]. we will always use p as the approx- 
imation to A in the harmonic projection method. 

Define r = (A — pl)y. Then (A — pl)y _L y. We replace r, the residual of standard Ritz 
pair, in ([7|) and by this new r and the standard Ritz vector y in ([8]) by the harmonic Ritz 
vector. We then solve the resulting (|T|) and ([8]) with the new r and y, respectively. With the 
solutions, we expand the subspaces analogous to Steps 4-5 of Algorithm 1. In such a way, 
we get the HSIRA and HJD methods. 

We should point out that a harmonic JD method had already been proposed as early as 
the standard one in [23] . where it is suggested to solve the following correction equation for 
u'j _L y' (suppose that y' H y ^ 0) 

with r = (A — pl)y and y' = (A — crl)y. It is different from (|8j) in the harmonic JD method 
described above, so its solution is different from that of ([8]) too. We prefer ([8]) since there is 
a potential danger that the oblique projector involved in (|lip may make it worse conditioned 
than (jHJ| considerably. 

In what follows we propose the RHSIRA and RHJD methods. We first compute the 
Rayleigh quotient p defined by (fTUj) and then seek a unit length vector y € V satisfying the 
optimality property 

\\{A-pI)y\\= min \\{A-pI)w\\ (12) 

ujGV, ||io||=1 

and use it to approximate x. So y is the best approximation to x from V with respect to 
p and the Euclidean norm. We call y a refined harmonic Ritz vector or more generally a 
refined approximate eigenvector since the pair (p, y) does not satisfy the harmonic projection 
([3]) any more. Since V forms an orthonormal basis of V, (112j) is equivalent to finding the unit 
length z € C m such that y = Vz with z satisfying 

11(^4 — pl)y\\ = min || (A — pI)Vz\\ 

z£C m , \\z\\ = l 

= \\{A-pI)Vz\\ 
= a mhl {(A- pI)V) . 
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We see that z is the right singular vector associated with smallest singular value of the matrix 
(A — pI)V, and y = Vz. So we can get z and \\(A — pI)Vz\\ by computing the singular value 
decomposition (SVD) of (A — pI)V. 
Similar to (llOj) . define 

p = y H Ay = z H H H z + a, (13) 

the Rayleigh quotient of A with respect to the refined harmonic Ritz vector y. Then the new 
residual r = {A — p)y _L y. Replace r and y by f and y in ([7]) and (|8j) and perform Steps 
3-5. Then we obtain the RHSIRA and RHJD methods, respectively. We will use (p, y) to 
approximate (A, x) in the methods as \\{A — pl)y\\ < \\(A — pl)y\\. 

Some important results have been established in |liyi2] for standard and refined projec- 
tion methods. The following two results in are directly adapted to the harmonic and 
refined harmonic projection methods: First, we have ||(^4 — pl)y\\ < \\(A — pl)y\\ unless the 
latter is zero, that is, the pair (p,y) is an exact eigenpair of A. Second, if \\(A — pl)y\\ ^ 
and there is another harmonic Ritz value close to p, then it may occur that 

\\(A-pI)y\\<^\\(A-pI)y\\. 

So y can be much more accurate than y as an approximation to the desired eigenvector x. 

For a general m-dimensional subspace V, two approaches are proposed in [9j[13] for com- 
puting z. Approach I is to directly form the cross-product matrix 

S = V H (A - P I) H (A - pI)V (14) 
= G+(F=p)H H + (a-p)H+\a-p\ 2 I, (15) 

which is Hermitian semi-positive definite. The desired z is the normalized eigenvector asso- 
ciated with the smallest eigenvalue of S, and ||(^4 — pl)y\\ is the square root of the smallest 
eigenvalue. Noticing that G and H are already available in the procedure of the harmonic 
projection, we can form S using at most 12m 2 flops by taking the worst case into account 
that all terms in f)15|) are complex. So, as a whole, we can compute z by the QR algorithm 
at cost of 0(m 3 ) + 12m 2 flops. 

Approach II is to first make the thin or compact QR decomposition (A — pI)V = QR and 
then make the SVD of the m x m triangular matrix R by 0(m 3 ) flops to get the smallest 
singular value and the associated right singular vector, which are \\(A — pl)y\\ and 5, respec- 
tively. If we use the Gram-Schmidt orthogonalization procedure with iterative refinement to 
compute the QR decomposition, then we will get a numerically orthonormal Q within a small 
multiple of the machine precision, which totally needs no more at most 4nm? flops generally 
if p is real. As a whole, 4nm 2 + 0(m 3 ) flops are needed. 

By comparison, we find that Approach I is, computationally, much more effective than 
Approach II. It can be justified [13] that Approach I is numerically stable for computing z 
provided that there is a considerable gap between the smallest singular value and the second 
smallest one of (^4 — pI)V . So in this paper, we use Approach I to compute y. 

Because of different right-hand sides, it is important to note that expanded subspaces are 
generally different for the SIRA, HSIRA and RSIRA methods whatever the linear systems are 
solved either exactly or approximately. The same is true for the JD, HJD and RHJD methods 
because of not only different right-hand sides but also different y in the correction equation 
([8]). This is different from shift-invert Arnoldi (SIA) type methods, i.e., the standard SIA, 
harmonic SIA, refined SIA and refined harmonic SIA, where the the updated subspaces are 
the same once the linear systems are solved exactly or approximately using the same inner 
solver with the same accuracy because the right-hand sides involved are always the currently 
newest basis vector at each outer iteration step. We will come back to this point in the end 
of Section [3] and show which method favors subspace expansion. 
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3 Accuracy requirements of inner iterations in HSIRA, HJD, 
RHSIRA and RHJD 



In this section, we review some important results in [T3] and apply them to establish the 
convergence theory of the inexact HSIRA, HJD and RHSIRA and RHJD. We prove that 
each inexact method mimics its exact counterpart well provided that all the inner linear 
systems are solved with only low or modest accuracy. We stress that in the presentation y is 
a general approximation to x. Returning to our methods, y is just the harmonic Ritz vector y 
in HSIRA and HJD and the refined harmonic Ritz vector y in RHSIRA and RHJD. Finally, 
we look into the way that each exact method expands the subspace and make a simple 
analysis, showing that HSIRA, HJD and RHSIRA, RHJD generally expand the subspace 
more effectively than SIRA and JD when computing interior eigenpairs, so that the former 
ones are more effective than the latter ones. This advantage conveys to the inexact HSIRA, 
HJD and RHSIRA, RHJD when each inexact method mimics its exact counterpart well. 
We can write the linear system ([7]) and the correction equation ([8]) as a unified form 

(A — al)u = a\y + a^A — al)y, (16) 

which leads to (J7J) in the HSIRA or RHSIRA method and © in the HJD or RHJD method 
when cti = a — u, a2 = 1 or a.\ = a — p, a.2 = 1 and a\ = — y H By , ^2 = 1 or ax = 
— y H By , OL2 = 1, respectively. Define B = (A — al)^ 1 . Then the exact solution u of (JX6|) is 

u = a\By + aiy ■ (17) 
Let u be an approximate solution of (|16p and its relative error be 

\\u — u\\ , 

(18) 



\\u\\ 

Then we can write 

u = u + e||«||/, (19) 

where / is the normalized error direction vector. Note that the direction of / depends on 
inner linear systems solves and accuracy requirements, but it shows no particular features 
in general. So sinZ(V, /) is generally moderate and not near zero. The vectors u and u 
are orthogonalized against V to get the (unnormalized) expansion vectors (J — Pv)u and 
(7 — Pv)u, respectively, where Py is the orthogonal projector onto V. Define the relative 
error 

_ ||(7 -P v )u- (7 -7VH 

11(1 ~PV)U\\ ' (20) 

which measures the relative difference between two expansion vectors (I—Py)u and (7— Py)u. 
The following result Theorem 3] establishes a compact bound for e in terms of e. 

Theorem 1. Define B = (A — crl)^ 1 and a = — ^ with u\ ^ 0, where a\ and «2 o^e given 
16]) . Suppose (-jr!— ,£) is a simple desired eigenpair of B and let (x,Xj_) be unitary. 



Then 



x H 



B [ x X A 



1 C H 

A — u 

L 



where c H = x H BX± and L = X^BX±. Assume that sinZ(V,/) ^ and a is not an 
eigenvalue of L and define 



sep(a,L) = |[(L-a7) _1 || _1 > 



7 



Then 

e < Ci, (21) 

where 

Q . = 2 W B W (22) 

' sep(a,L)sinZ(V,/) V ' 

with a = or a = in SIRA type methods and a = y H By or y H By in JD type methods. 

If a is a fairly approximation to and y^— is not close to any eigenvalue of L, then 
sep(a, L) is not small and is actually of 0(||5||) provided that the eigensystem of B is not ill 
conditioned. In this case, C is moderate as sinZ(V, /) is moderate, as commented previously. 
For a given e, we should note that the bigger sinZ(V, /) is, the bigger e is. So it is a lucky 
event if sinZ(V, /) is big as it means that we need to solve inner linear systems with less 
accuracy e. 

Below we discuss how to determine e such that the inexact methods mimic their exact 
counterparts very well, that is, each inexact method and its exact counterpart use comparable 
or almost the same outer iterations to achieve the convergence. The following important result 
is proved in |14j . which forms the basis of our analysis. 

Theorem 2. Define 

V+ = V©span{(7-P V M, V+ = V span {(/ - P v )u} , 
where u and u are the exact and inexact solutions of the linear system (|16fl . respectively, and 

sinZ(V + ,x) ~ smZ(V + ,x) 

smZ{V,x) ' sinZ(V,x) ' [ ' 



Then we have 



Furthermore, if r < 1, then 



smZ(V+,x) 8 . , 

■ / v = 7 € 1-t,1 + t. 25 
smZ(l/+,x) d 

(]23p measures one step subspace improvements of the exact and inexact subspace ex- 
pansions, respectively. ([25]) indicates that to make sinZ(V + ,x) sinZ(V+,x), r should be 
small. Note that the difference of the upper and lower bounds is 2r. So a very small r only 
improves the bounds very marginally, and a fairly small r, e.g., r £ [0.001,0.01], is enough 

since we have G [0.999, 1.001] or [0.99, 1.01] and the lower and upper bounds differ 

marginally, which means that the two subspaces V+ and V+ are of comparable or almost 
the same quality for a fairly small r when computing (A, x). Precisely, with respect to the 
two subspaces the approximations of the desired (X,x) obtained by the same type, i.e., the 
Rayleigh-Ritz method, its harmonic variant or the refined harmonic variant, should generally 
have the comparable or almost the same accuracy. By (|24p . we have 

T X 5 , , 

e=— (26) 

As it has been proved in [14] , the size of 5 crucially depends on the eigenvalue distribution 
of B. Generally speaking, the better is separated from the other eigenvalues of B, the 
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smaller 5 is. Conversely, if is poorly separated from the others, 5 may be near to 
one; see [13] for a descriptive analysis. 5 is an a-priori quantity and is unknown during 
computation. But generally we should not expect that a practical problem is very well 
conditioned, that is, 5 is not much less than one; otherwise, the exact methods will generally 
find (A, x) by using only a very few outer iterations. For example, suppose that the initial 
V is one dimensional and 5 = 0.1 at each outer iteration. Then the updated sinZ(V+,x) 
is no more than 10 -10 after ten outer iterations and V+ is a very accurate subspace for 
computing x. So in practice, we assume that 5 is not much less than one, say no smaller than 
0.2. As we have argued, a fairly small r € [0.001,0.01] should make the inexact methods 
behave very like their exact counterparts. As a result, by (|26p and the argument that a fairly 
small r 6 [0.001,0.01] is generally enough , to make the inexact methods mimic their exact 
counterparts, it is reasonable to take 



which is also suggested in [14] for the standard SIRA and JD methods. For r € [0.001,0.01], 
if e defined by (|27p is unfortunately bigger than that in (|26p . the the inexact methods may 
use more outer iterations and may not behave very like the exact methods. Even so, however, 
since solving inner linear systems with e in (|27p is cheaper than with considerably smaller e, 
the overall efficiency of the inexact methods may still be improved. 

Finally, for the exact methods, let us, qualitatively and concisely, show which updated 
subspace V+ is better and which method behaves more favorable. Keep (|1T|) and y G V in 
mind. The (unnormalized) expansion vector is 



So we actually use add By to V to get the expanded subspace V+. It is easily justified 
that the better y approximates x, the better does By in direction. So for a more accurate 
approximate eigenvector y, By and thus the expanded subspace V+ contain richer information 
on the eigenevector x. As we have argued, since the harmonic Ritz vector is a more reliable 
and regular approximation to x while the Ritz vector may not be for the interior eigenvalue 
problem, so the exact HSIRA and HJD may expand the subspaces better and more regularly 
than the exact SIRA and JD. Furthermore, since the refined harmonic Ritz vector is generally 
more and can be much more accurate than the harmonic Ritz vector, the exact RHSIRA and 
RHJD, in turn, generate better subspaces than the exact HSIRA and HJD at each outer 
iteration. As a consequence, HSIRA and HJD are expected to converge more regularly and 
use fewer outer iterations than SIRA and JD while RHSIRA and RHJD may use the fewest 
outer iterations among all the exact methods for the interior eigenvalue problem. For the 
inexact methods, provided that the selection of e makes each inexact method mimic its exact 
counterpart well, the inexact HSIRA, HJD and RHSIRA, RHJD are advantageous to the 
inexact SIRA and JD. Numerical experiments will confirm our expectations. 

4 Restarted algorithms and stopping criteria for inner solves 

Due to the storage requirement and computational cost, it is generally necessary to restart 
the methods to avoid large steps of outer iterations. We describe the restarted HSIRA/HJD 
algorithms as Algorithm [2] and their refined variants as Algorithm [3j 

We consider some practical issues on two algorithms. For outer iteration steps m = 
1,2,... ,M max during the current cycle, suppose {p^ m \ j/ m ') is used to approximate the de- 
sired eigenpair (A, x) of A at the m-th outer iteration. As done commonly in the literature, 



£ G [10~ 4 , 10~ 3 ] 



(27) 



(/ - P v )u 



{I - P v )By. 
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Algorithm 2 Restarted HSIRA/HJD algorithm with the fixed target a 

Given the target a and a user-prescribed convergence tolerance tol, suppose the columns 
of V form an orthonormal basis of an initial subspace V and let M max be the maximum of 
outer iterations allowed, 
repeat 

1. Compute (update) H = V H (A - aI) H V and G = V H {A - aI) H (A - aI)V. 

2. Let (/i, z) be an eigenpair of the matrix pencil (H, G), where fi = A. 

3. Compute the Rayleigh quotient p = z H H H z + a and the harmonic Ritz vector 
y = Vz. 

4. Compute the residual r = Ay — py. 

5. In HSIRA or HJD, solve the inner linear system 

(A-aI)u = r or (I - yy H )(A - al)(l - yy H )u = -r. 

6. Orthonormalize u against V to get the expansion vector v. 

7. If dim(V) < M max , set V = [ V v ] ; otherwise, V = V new . 
until llrll < tol 



Algorithm 3 Restarted RHSIRA/RHJD algorithm with the fixed target a 

Given the target a and a user-prescribed convergence tolerance tol, suppose the columns 
of V form an orthonormal basis of an initial subspace V and let M max be the maximum of 
outer iterations allowed, 
repeat 

1. Compute (update) H = V H {A - aI) H V and G = V H (A - aI) H (A - aI)V. 

2. Let (fi, z) be an eigenpair of the matrix pencil (H, G), where fi = A. 

3. Compute the Rayleigh quotient p = z H H H z + a. 

4. Form the cross-product matrix S = G + (a — p)H H + (a — p)H + \o~ — p\ 2 I and 
compute the eigenvector z of S associated with its smallest eigenvalue. 

5. Compute the new Rayleigh quotient p = z H H z, + a and the refined harmonic 
Ritz vector y = Vz. 

6. Compute the residual r = Ay — py. 

7. In RHSIRA or RHJD, solve the inner linear system 

(A - al)u = r or {I - yy H )(A - al)(l - yy H )u = -r. 

8. Orthonormalize u against V to get the expansion vector v. 

9. If dim(V) < M max , set V = [ V v ]; otherwise, V = V new . 
until llrll < tol 



we simply take the new starting vector 

\r — ,.(M max ) 
''new — y 

to restart the algorithms. In practice, if y( Mmax ) is complex conjugate, we take their real and 
imaginary parts, normalize and orthonormalize them to get an orthonormal V ncw of column 
two. We mention that a thick restarting technique [24] may be used, in which, besides y( Mmax ) ) 
one retains the approximate eigenvectors associated with a few other approximate eigenvalues 
closest to p( Mmax ) and orthonormalizes them to obtain V^ cw . We will not report numerical 
results with thick restarting. 
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Note that all the methods under consideration do not have residual monotonically de- 
creasing properties as outer iteration steps increase. Therefore, it may be possible for both 
Algorithm [2] and Algorithm [3] to take bad restarting vectors. This is indeed the case for 
the standard Rayleigh-Ritz method when computing interior eigenpairs, as has been widely 
realized, e.g., |26j. However, as was stated previously, the harmonic and refined harmonic 
methods are more reliable to select correct and good approximations to the desired eigen- 
pairs. So they are more suitable and reliable for restarting than the standard Ritz vector for 
the interior eigenvalue problem. As a consequence, Algorithm [2] and especially Algorithm [3] 
converge more regularly than the restarted standard SIRA and JD. Our numerical examples 
will confirm this and illustrate that the refined harmonic Ritz vectors are better than the 
harmonic Ritz vectors for restarting and Algorithm [3] converges faster than Algorithm [2j 

Our ultimate goal is to get a practical estimate for e for a given e. The a-priori bound (|2ip 
forms the basis for this. Let v\, U2, ■ ■ ■ , v m are the harmonic Ritz values, i.e., the eigenvalues 
of the pencil (H m , G m ) at the m-th outer iteration during the current cycle (here we add the 
subscript m to H and G in the algorithms) and assume that v\ is used to approximate the 
desired eigenvalue A. We simply estimate ~ \ P -a\ > wnere P ls the Rayleigh quotient of 
A with respect to the harmonic Ritz vector in HSIRA and HJD and the refined harmonic 
Ritz vector in RHSIRA and RHJD. For sep(ct, L), we replace a, an approximation to jz^, 
by -^-^ accordingly and then estimate 

i i 

sep(a, L) min 



2,3,...,m p — a — a 

Finally, replace sinZ(V,/) by its maximum one. Combining all these together, we get the 
following estimate C of C in (|2ip : 

— , m > 1 

-P , (28) 

1, m = 1 



2 max 
C' = < i=2,3,...,i 



This is analogous to what is done in [14], and the difference is that we here use the harmonic 
Ritz values z^, . . . , v m to replace the standard Ritz values used in p3] as approximations to 
some eigenvalues of A other than A. Denote by eg and ej practical e's used in the HSIRA, 
RHSIRA and HJD, RHJD algorithms, respectively. Recall that bound (|2ip is compact. Then 
we take 

es = ej = e = Ce. (29) 

For a fairly small e, we may have e > 1 in case C is big, which will make u no sense as an 
approximation to u. In order to make u have some accuracy, we propose using 

e = min{C"e,0.l} , (30) 

which shows that e is comparable to e in size whenever i is small and C is moderate. 



5 Numerical experiments 

All the numerical experiments were performed on an Intel (R) Core (TM)2 Quad CPU 
Q9400 2.66GHz with the main memory 2 GB using Matlab 7.8.0 with the machine precision 
£mach = 2.22 x 10~ 16 under the Linux operating system. All the test examples are difficult 
in the sense that the desired eigenvalue is clustered with some other eigenvalues of A for the 
given target a. We aim to show four points. First, regarding restarts of outer iterations, 
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for fairly small e = 10" 3 and 1(T 4 , the restarted inexact HSIRA/HJD and RHSIRA/RHJD 
algorithms behave (very) like the exact counterparts. Second, regarding outer iterations, 
Algorithms [2H3] are much more efficient than the restarted standard SIRA and JD for the 
same e's, and Algorithm [3] is the best. Third, regarding total inner iterations and overall 
efficiency, Algorithm [3] is considerably more efficient than Algorithm [2] and the restarted 
inexact standard SIRA and JD perform very poorly and often fail to converge. Fourth, each 
SIRA type algorithm is equally as efficient as the corresponding JD type algorithm for the 
same e. 

At the m-th outer iteration step of the inexact HSIRA or HJD method, we have H m = 

( m ) ( m )\ ■ 1 1 

, m 



V*(A - aI) H V m and G m = V£ (A - aI) H (A - aI)V m . Let {y\ n \ V^), i = 1,2 
be the harmonic 
in mind. We use 



be the harmonic Ritz pairs, labeled as \dy — cr| < |i>2 ~ a \ 5; ' ' ' < \ u m^ ~ a \- Keep (fit) 



(p M y M) = ^Y H ^z[ m) + a, V m z[ m) ) (31) 

to approximate the desired eigenpair (A, x) of A. If the RSIRA and RHJD methods are used, 
we form the cross-product matrix S m by f)15|) and compute its eigenvector associated with 
the smallest eigenvalue. We then compute the refined harmonic Ritz vector y( m ' = V m z^ 
and the Rayleigh quotient p( m > defined by (fT3|) . Let r m = Ay( m ^ — p( m )y( m ) . We stop the 
algorithms if 

||r m || < tol = max{p||i,l} x 10~ 12 . (32) 
In the algorithms, we stop inner iterations at outer iteration m when 
\\r m -(A- al)u\\ || - r m - (I - y ( - m \y^) H )(A - al)(l - y( m )(y( m )) H )n|| 



< e, (33) 

II ' m || || ' m || 

where e is defined by ([3D]) for a given e. We will denote by SIRA(e), JD(e), HSIRA (e), 
HJD(e) and RHSIRA(e), RHJD(e) the inexact algorithms with the given e. 

In the "exact" SIRA and JD type algorithms, we stop inner iterations when (|33p is satisfied 
with e = 10" 14 . 

All the test problems are from Matrix Market pQ. For each inner linear system, we 
always took the zero vector as an initial approximate solution and solved it by the right- 
preconditioned restarted GMRES(30) algorithm. Each algorithm starts with the normalized 
vector 7^(1; 1) • • • i ty H ■ At the m-th outer iteration of current restart, m = 1,2,..., M max , 
for the correction equations in the JD type algorithms we used the preconditioner 

M m = (l - y {m) (y {m) ) H ) M(l- y {m) (y {m) ) H ) (34) 

suggested in [26] , where M ~ A — a I is the incomplete LU preconditioner for the SIRA type 
algorithms. For all the algorithms, the maximum steps of outer iterations are 30 per restart. 
An algorithm signals the failure if it did not converge within Max restarts. 

In the experiments, each figure consists of three subfigures, in which the top subfigure 
denotes outer residual norms versus outer iterations of the first cycle before restart, the 
bottom-left subfigure denotes outer residual norms versus restarts and the bottom-right sub- 
figure denotes the numbers of inner iterations versus restarts. 

The top subfigures exhibit the qualities of the standard, harmonic and refined harmonic 
Ritz vectors for restarting. The bottom two subfigures depict the convergence processes of 
the exact SIRA type algorithms and their inexact counterparts with e = 10 -3 , showing the 
convergence behavior of the algorithms and the local efficiency per restart, respectively. 
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In all the table, denote by Irestart the number of restarts to achieve the convergence, by 
hnner the total number of inner iterations and by Prj.i the percentage of the times /o.i that 
e = 0.1 is used in the total number of outer iterations Iouter, i-e., 

P - /qi 
-nu — 7 • 



1 outer 



Note that Ii nner equals the total products of A and vectors in the restarted GMRES algorithm. 
It is a reasonable measure of the overall efficiency of all the algorithms used in the experiments. 



Example 1. This unsymmetric eigenproblem M80PI_n of n = 4182 arises from real power 
system models pQ. We test this example with a = 0.05 + 0.5i and Max = 500. The computed 
eigenvalue is A ~ —6.9391 x 10~ 5 + 5.0062 x 10 _1 i. The preconditioner M is the incomplete 
LU factorization of A — al with drop tolerance 10" 1 . Figure[T]and Table Q] display the results. 




5 10 15 20 25 30 

outer iterations 




Figure 1: Example 1. M80PLn with a = 0.05 + 0.5i. 



We see from the top subfigure of Figure[T]that for e = 10~ 3 the inexact SIRA, HSIRA and 
RHSIRA behaved very like their corresponding exact counterparts in the first cycle but the 
exact and inexact RHSIRA are more effective than the exact and inexact HSIRA. It clearly 
shows that the exact SIRA and SIRA(10 -3 ) are the poorest and considerably poorer than 
the corresponding HSIRA and RHSIRA. There are two reasons for this. The first is that 
harmonic and refined harmonic Ritz vectors are more reliable than standard Ritz vectors 
for computing interior eigenvectors. The second is that the harmonic and refined harmonic 
Ritz vectors favor subspace expansions. We, therefore, expect that the restarted HSIRA 
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Accuracy 


Algorithm 


-^restart 


firmer 


Po.i 


Algorithm 


1 restart 


dinner 


Po.i 




SIRA 


Max 


32142 


0% 


JD 


Max 


33823 


0% 


e = 1CT 3 


HSIRA 


383 


35032 


0% 


HJD 


371 


33988 


0% 




RHSIRA 


79 


4981 


89% 


RHJD 


79 


4958 


89% 




SIRA 


Max 


47166 


0% 


JD 


Max 


47785 


0% 


£ = 10- 4 


HSIRA 


376 


46470 


0% 


HJD 


374 


46231 


0% 




RHSIRA 


74 


7032 


0% 


RHJD 


75 


7148 


0% 




SIRA 


Max 


217941 




JD 


Max 


217936 




exact 


HSIRA 


368 


171183 




HJD 


368 


171182 






RHSIRA 


76 


35515 




RHJD 


76 


35513 





Table 1: Example 1. M80PLn with a = 0.05 + 0.5i. 



and RHSIRA can be much more efficient and converge much faster than the restarted SIRA 
when computing interior eigenpairs since we use more reliable and possibly accurate harmonic 
and refined harmonic Ritz vectors as restarting vectors. Furthermore, as far as restarts are 
concerned, the restarted RHSIRA outperforms the restarted HSIRA very substantially. We 
observed very similar convergence behavior for the JD type algorithms and thus have the same 
expectations on the restarted JD type algorithms. These expectations are indeed confirmed 
by numerical experiments, as shown by I res tart' s i n the table and figure. 

We explain the table and figure in more details. From the bottom-left subfigure, regarding 
outer iterations, we see that both the exact and inexact restarted SIRA did not converge 
within 500 restarts while HSIRA and RHSIRA worked well and the inexact methods behaved 
very like their corresponding exact ones. The restarted RHSIRA and RHJD were the fastest 
and five times as fast as the restarted HSIRA and HJD, respectively. As is expected, the 
table confirms that, for the same e, each SIRA type algorithm and the corresponding JD 
type algorithm behaved very similar and were almost indistinguishable. 

Regarding the overall efficiency, the exact HSIRA, HJD, RHSIRA and RHJD each used 
about 460 inner iterations per restart. In contrast, HSIRA(10" 3 ) and HSIRA(10" 4 ) used 
almost constant inner iterations each restart, which were about 100 and 110 per restart, 
respectively, and RHSIRA(10- 3 ) and RHSIRA(10~ 4 ) used roughly 60 and 100 inner itera- 
tions per restart, respectively. The same observations are true for the JD type algorithms. 
These experiments demonstrate that modest e = 10 -3 , 10 are enough to make the in- 
exact restarted algorithms mimic their exact counterparts. We find that the SIRA(10~ 4 ) 
and JD(10~ 4 ) type algorithms used almost the same outer iterations as the SIRA(10 -3 ) and 
JD(10~ 3 ) type ones, but the latter consumed considerably fewer total inner iterations and 
improved the overall efficiency substantially. So, smaller e's are not necessary for this example 
as they cannot reduce outer iterations further and may cost much more inner iterations. 

In addition, we see from Table [T] that 89 percent of inner linear systems in RHSIRA(10 -3 ) 
were actually solved with the accuracy requirement e = 0.1. This means that though most 
of C"s in (|28|) are big, it suffices to solve all the inner linear systems with low accuracy 0.1. 

Example 2. This unsymmetric eigenproblem M80PI_n of n = 4182 arises from real power 
system models [T]. We test test this example with a = 0.4 + 1.3i and Max = 1000. The 
computed eigenvalue is A ~ —2. 8266 x 10 _4 -|-1.3068i. The preconditioner M is the incomplete 
LU factorization of A — al with drop tolerance 10 _1 . Figure[2]and Table [2] report the results. 

The top subfigure of Figure [2] is similar to that of Figure [H showing that harmonic 
and especially refined harmonic Ritz vectors are more suitable for expanding subspaces and 
restarting for the interior eigenvalue problem. The difference is that this problem is more 
difficult than Example 1 since, for the 30-dimensional subspace in the first cycle, the residual 
norms decrease more slowly and approximate eigenpairs are less accurate than those for 
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Figure 2: Example 2. S80PLn with a = 0.4 + 1.3i. 



Accuracy 


Algorithm 


-^restart 


dinner 


Pq.i 


Algorithm 


^restart 


dinner 


Po.i 




SIRA 


Max 


43020 


6% 


JD 


Max 


43279 


5% 


e = 1CT 3 


HSIRA 


Max 


67645 


0% 


HJD 


Max 


67686 


0% 




RHSIRA 


926 


29774 


96% 


RHJD 


946 


30392 


96% 




SIRA 


Max 


82712 


0% 


JD 


Max 


82335 


0% 


i = icr 4 


HSIRA 


Max 


104141 


0% 


HJD 


Max 


104102 


0% 




RHSIRA 


945 


31323 


94% 


RHJD 


937 


31068 


93% 




SIRA 


Max 


377424 




JD 


Max 


377426 




exact 


HSIRA 


Max 


382069 




HJD 


Max 


381359 






RHSIRA 


699 


271142 




RHJD 


684 


263865 





Table 2: Example 2. S80PLn with a = 0.4 + 1.3i. 



Example 1. So it is expected that the restarted algorithms converge more slowly and use 
more outer iterations Irestart® than for Example 1. We observe from the figure that the 
convergence curves of the three exact algorithms SIRA, HSIRA and RHSIRA essentially 
coincide with those of their inexact variants with e = 10~ 3 . So smaller e doe not help, 
rather it makes the algorithms may waste much more inner iterations. We observed similar 
convergence phenomena for the JD type algorithms. 

Tableland Figure [2] shows that the restarted exact and inexact SIRA, JD, HSIRA and 
HJD all failed to converge within 1000 restarts. Furthermore, it is deduced from Figure[2]that 
the restarted exact and inexact SIRA and JD appear impossible to converge at all since their 
convergence curves were irregular, oscillated and had no decreasing tendency. The restarted 
exact and inexact SIRA and JD behaved regular but converged too slowly. In contrast, 
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the restarted exact and inexact RHSIRA and RHJD with e = 10 -3 and 10 -4 converged. 
Furthermore, the restarted RHSIRA and RHJD with i = 10~ 3 and e = 10 -4 behaved very 
similar and used almost the same outer iterations. Both of them mimic the restarted exact 
RHSIRA and RHJD fairly well. So, smaller e cannot reduce outer iterations substantially. 

Regarding the overall efficiency, the inexact restarted algorithms improved the situation 
tremendously. As in Example 1, the exact RHSIRA and RHJD still needed about 370 inner 
iterations per restart, much more than 40 inner iterations that were used by the inexact 
RHSIRA and RHJD for e = 10~ 3 and 10 -4 per restart. As a whole, Imner's illustrate 
that the inexact algorithms with these two e's were about eight times faster than the exact 
algorithms, a striking improvement. Besides, note that the number of inner linear systems 
that were solved with lower accuracy e = 0.1 in RHSIRA(10~ 3 ) were 96%, while those in 
HSIRA(10~ 3 ) and SIRA(10~ 3 ) were 0% and 6%, respectively. This is why RHSIRA(10~ 3 ) 
used fewer inner iterations than SIRA(10" 3 ) and HSIRA(10~ 3 ) per restart, as shown in the 
bottom-right subfigure. 

Although the restarted exact SIRA and HSIRA both failed for this example, the reasons 
may be completely different. In the top subfigure, we find that the convergence curves of 
SIRA bulged at the last few steps (25 ~ 30). This is harmful to restarting as it is very possible 
to get an unsatisfying restarting vector once the method bulged at the very last step. In the 
bottomdeft subfigure, it is seen that the convergence curves of restarted exact and inexact 
SIRA were irregular while the convergence curves of HSIRA decreased smoothly though the 
algorithm converged quite slow. So we can infer that HSIRA did not take bad restarting 
vectors. Remarkably, the figure and table tell us that RHSIRA converged much faster. This 
is due to the fact that the refined harmonic Ritz vector y can be much more accurate than 
the corresponding harmonic Ritz vector y, so that restarting vectors were better and the 
subspaces generated were more accurate as restarts proceeded. 

Example 3. This unsymmetric eigenproblem dw4096 of n = 8192 arises from dielectric 
channel waveguide problems [I]. We test this example with a = — 24 and Max = 500. 
The computed eigenvalue is A ~ —30.217. The preconditioner M is the incomplete LU 
factorization of A — al with drop tolerance 10~ 3 . Figure [3] and Table [3] display the results. 



Accuracy 


Algorithm 


^restart 


dinner 


Po.i 


Algorithm 


^restart 


dinner 


Po.i 




SIRA 


Max 


15722 


45% 


JD 


Max 


17118 


45% 


e = 1CT 3 


HSIRA 


306 


9473 


85% 


HJD 


300 


9246 


85% 




RHSIRA 


168 


5019 


96% 


RHJD 


151 


4540 


96% 




SIRA 


Max 


20586 


14% 


JD 


Max 


20794 


14% 


e = 1CT 4 


HSIRA 


289 


9251 


49% 


HJD 


282 


9091 


49% 




RHSIRA 


156 


4722 


94% 


RHJD 


197 


5952 


95% 




SIRA 


Max 


87410 




JD 


Max 


89435 




exact 


HSIRA 


255 


44906 




HJD 


259 


45480 






RHSIRA 


127 


23351 




RHJD 


139 


25419 





Table 3: Example 3. dw4096 with a = —24. 



Figure [3] indicates that the exact and inexact SIRA oscillated sharply in the first cycle, the 
exact and inexact HSIRA improved the situation very significantly but still did not behave 
quite regularly, while the exact and inexact RHSIRA converged very smoothly after first a 
few outer iterations and considerably faster than the HSIRA. Poor Ritz vectors further led to 
poor subspace expansion vectors, generating a sequence of poor subspaces. We observed very 
similar phenomena in the JD type algorithms. This means that the standard Ritz vectors 
were not suitable for restarting, and the harmonic Ritz vectors were much better but inferior 
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Figure 3: Example 3. dw4096 with a = —24. 



to the refined harmonic Ritz vectors for restarting. So it is expected that the restarted 
exact and inexact SIRA and JD algorithms may not work well, but the restarted HSIRA 
and HJD may work much better than the former ones and the restarted RHSIRA and RHJD 
outperform the HSIRA and HJD very considerably. 

The above expectations are confirmed by and Table [3] and the bottom-left subfigure of 
Figure [3l It is seen that the restarted exact and inexact SIRA and JD algorithms failed to 
converge within 500 restarts while HSIRA and HJD solved the problem successfully and the 
refined RHSIRA and RHJD were twice as efficient as the harmonic algorithms, as indicated 
by IrestartS- Furthermore, we see the inexact algorithms with e = 10~ 3 , 10 -4 exhibited very 
similar convergence behavior and mimic the exact algorithms very well. This confirms our 
theory that modest e is generally enough to make the inexact SIRA and JD type algorithms 
behave very like the corresponding exact counterparts and smaller e is not necessary. We 
also find that each exact SIRA type algorithm is as efficient as the corresponding JD type 
algorithm for the same e. 

Regarding the overall performance, the comments on Examples 1-2 apply here analo- 
gously. For E = 10" 3 and 10~ 4 , the inexact HSIRA, RHSIRA, HJD and RHJD were four 
times faster than the corresponding exact algorithms. Furthermore, the restarted exact and 
inexact RHSIRA and RHJD were twice as fast as the corresponding HSIRA and HJD algo- 
rithms for the same e. 

Example 4. This unsymmetric eigenproblem dw8192 of n = 8192 arises from dielectric 
channel waveguide problems [TJ. We test this example with a = —60 and Max = 1000. 
The computed eigenvalue is A ~ —87.795. The preconditioner M is the incomplete LU 
factorization of A — al with drop tolerance 10 -3 . Figure U] and Tables U] report the results. 
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Figure 4: Example 4- dw8192 with a = —60. 
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^restart 
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Po.i 


Algorithm 


^restart 


^inner 


Po.i 




SIRA 


Max 


29262 


39% 


JD 


Max 


32222 


39% 


e= lCT 3 


HSIRA 


542 


15953 


76% 


HJD 


402 


11848 


81% 




RHSIRA 


372 


10813 


96% 


RHJD 


317 


9221 


96% 




SIRA 


Max 


38502 


12% 


JD 


Max 


39413 


12% 


e = 1(T 4 


HSIRA 


555 


17416 


51% 


HJD 


497 


15691 


57% 




RHSIRA 


386 


11595 


93% 


RHJD 


429 


12884 


93% 




SIRA 


Max 


145435 




JD 


Max 


151453 




exact 


HSIRA 


615 


89636 




HJD 


617 


90056 






RHSIRA 


387 


59170 




RHJD 


383 


59166 





Table 4: Example 4. dw8192 with a = -60. 



All the algorithms gave results that are typically similar to those corresponding counter- 
parts obtained for Example 3. The comments and conclusions are very analogous, and we 
thus omit them. 



6 Conclusions 

The standard SIRA and JD methods may not work well for the interior eigenvalue problem. 
The standard Ritz vectors may converge irregularly, and it is also hard to select correct Ritz 
pairs to approximate the desired interior eigenpairs. So the Ritz vectors may not be good 
choices for restarting, causing that restarted algorithms may perform poorly. Meanwhile, 
the Ritz vectors appear to expand subspaces poorly. In contrast, the harmonic Ritz vectors 
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are more regular and more reliable approximations to the desired eigenvectors, so that they 
may expand the subspaces better and generate more accurate a sequence of subspaces when 
restarting the methods. Due to the optimality, the refined harmonic Ritz vectors are generally 
more accurate than the harmonic ones and are better choices for expanding subspaces or 
restarting the methods. Most importantly, we have proved that the harmonic and refined 
SIRA and JD methods generally mimic their exact counterparts well provided that all inner 
linear systems are solved with only low or modest accuracy. To be practical, we have presented 
the restarted harmonic and refined harmonic SIRA and JD algorithms. Meanwhile, we have 
designed practical stopping criteria for inner iterations. 

Numerical experiments have confirmed our theory. They have indicated that the restarted 
harmonic and refined harmonic SIRA and JD algorithms are much more efficient than the 
restarted standard SIRA and JD algorithms for the interior eigenvalue problem. Furthermore, 
the refined harmonic algorithms are much better than the harmonic ones. Importantly, the 
results have demonstrated that each inexact algorithm behaves like its exact counterpart 
when all inner linear systems are solved with only low or modest accuracy. In addition, 
numerical results have confirmed that each exact or inexact SIRA type algorithm is equally 
as efficient as the corresponding JD type algorithm for the same e. 

Our algorithms are designed to compute only one interior eigenvalue and its associated 
eigenvector. This is a special nature of SIRA and JD type methods: they compute only 
one eigenpair each time. If more than one eigenpairs are of interest, we may extend the 
algorithms to this situation in a number of ways. For example, we may introduce some 
deflation techniques [251126] into them. Also, we may change the target a to a new one, to 
which the second desired eigenvalue is the closest, and apply the algorithms to find it and 
its associated eigenvector. Proceed this way until all the desired eigenpairs are found. Such 
kind of algorithms is under developments. 
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